setwd("/Users/davidcarter/Dropbox/Systemic Crisis Rep/")


library(readstata13);library(foreign);library(lme4)


data<-read.dta13("systemic_data_compiled.dta") 




## Model 1 (Produce Figure 4, first plot) ##


lm1<-lmer(make_claim_onset_f~(sd_density-1|year)+mean_density+shape_area_10+ x+ y+onsetyrs+onsetyrs2+onsetyrs3+as.factor(dyad_m),data)


hi_ft1<-ranef(lm1)$year[,1]+ 1.96*sqrt(VarCorr(lm1)[[1]][1,1])
lo_ft1<-ranef(lm1)$year[,1]- 1.96*sqrt(VarCorr(lm1)[[1]][1,1])


yrs<-seq(1816,2002)


plot(1816:2002,ranef(lm1)$year[,1],pch=16,col="light grey",axes=F,xlab="",ylim= c(min(lo_ft1),max(hi_ft1)),ylab=expression(beta[t]),main="Model 1" )



points(seq(1816,2002)[lo_ft1 > 0], ranef(lm1)$year[lo_ft1 > 0,1], pch=16,col="dark red")

for(i in 1:length(yrs)){
  
  
  lines(c(yrs[i],yrs[i]),c(hi_ft1[i],lo_ft1[i]) )
  
  
  if(lo_ft1[i] > 0){
          lines(c(yrs[i],yrs[i]),c(hi_ft1[i],lo_ft1[i]), col="dark red",lwd=2 )
    
    }
  
  
}


axis(2)

axis(1,at=c(1816,1848,1866,1914,1938,1991,2012))





## Model 2 (Produce Figure 4, second plot) ##



lm2<-lmer(make_claim_onset_f~(sd_density-1|year)+mean_density+major1 +major2 +defense +both_democracy+shape_area_10+ x+ y+onsetyrs+onsetyrs2+onsetyrs3+as.factor(dyad_m),data)


hi_ft2<-ranef(lm2)$year[,1]+ 1.96*sqrt(VarCorr(lm2)[[1]][1,1])
lo_ft2<-ranef(lm2)$year[,1]- 1.96*sqrt(VarCorr(lm2)[[1]][1,1])


plot(1816:2002,ranef(lm2)$year[,1])





plot(1816:2002,ranef(lm2)$year[,1],pch=16,col="light grey",axes=F,xlab="",ylim= c(min(lo_ft1),max(hi_ft1)),ylab=expression(beta[t]),main="Model 2" )



points(seq(1816,2002)[lo_ft2 > 0], ranef(lm2)$year[lo_ft2 > 0,1], pch=16,col="dark red")

for(i in 1:length(yrs)){
  
  
  lines(c(yrs[i],yrs[i]),c(hi_ft2[i],lo_ft2[i]) )
  
  
  if(lo_ft2[i] > 0){
    lines(c(yrs[i],yrs[i]),c(hi_ft2[i],lo_ft2[i]), col="dark red",lwd=2 )
    
  }
  
  
}


axis(2)

axis(1,at=c(1816,1848,1866,1914,1938,1991,2012))





